# fundraising paper -- supplemental analysis estimating individual level 
# marginal component effects using hierarchical model

library(tidyverse)
library(brms)
library(tidybayes)
library(tidyr)
library(stringr)
library(kableExtra)
library(ggthemes)

overleaf.dir <- "~/Dropbox/Apps/Overleaf/Stepping Up the Political Ladder Appendix/"
  
  
conjoint <- readRDS("Modified Data/Clean_reshaped.RDS")


## make text bigger!
theme_set(theme_bw() + 
            theme(panel.spacing = unit(0, "in"), 
                  panel.grid = element_blank(),
                  text = element_text(size = 14))
)


# outcome variable- interest in running, rescaled so equally spaced from 0 to 1 
conjoint = conjoint %>% 
  mutate(interest_bin = as.numeric(interest > 0),
         interest_num = interest - min(interest,na.rm=TRUE)) %>% 
  mutate(interest_num = interest_num / max(interest_num,na.rm=TRUE))



# fit models ---------------------------------------------------------------

# estimates random slopes by respondent, allowing for individual-level heterogeneity. 
# see Jensen, Marble, Scheve, Stasavage (2021, PSRM) for explanation

interest_reg_hier = brm(interest_num ~ 
                        salary + fundraising + advertising + opp_exp + opp_ideol + 
                        (1 + salary + fundraising + advertising + opp_exp + opp_ideol || responseid), data = conjoint,
                       backend = "cmdstanr", iter = 1000, cores = 4, chains = 4,
                       control = list(max_treedepth = 13, adapt_delta = .95))



winprob_reg_hier = brm(winprob ~ 
                          salary + fundraising + advertising + opp_exp + opp_ideol + 
                          (1 + salary + fundraising + advertising + opp_exp + opp_ideol || responseid), data = conjoint,
                        backend = "cmdstanr", iter = 1000, cores = 4, chains = 4,
                        control = list(max_treedepth = 13, adapt_delta = .95))






# Results for interest outcome --------------------------------------------


# extract point estimates of individual-level effects
imces_pes <- coef(interest_reg_hier)[[1]][,1,] # point estimates of IMCEs
rids <- rownames(imces_pes) 
imces_pes <- as_tibble(imces_pes) %>% 
  mutate(responseid = rids) %>% 
  pivot_longer(-responseid) 

# tidy data frame for plotting  
imces_pes <- imces_pes %>% 
  mutate(
    attribute = case_when(
      grepl("salary", name) ~ "Salary",
      grepl("fundraising", name) ~ "Fundraising",
      grepl("opp_exp", name) ~ "Opponent's Experience",
      grepl("opp_ideol", name) ~ "Opponent's Ideology",
      grepl("advertising", name) ~ "Advertising"),
    variable = case_when(
      name == "Intercept" ~ "Intercept",
      name == "salary$50000" ~ "$50,000",
      name == "salary$80000" ~ "$80,000",
      name == "fundraising$100000" ~ "$100,000",        
      name == "fundraising$300000" ~ "$300,000",
      name == "advertisingMostlynegativeads" ~ "Mostly negative ads",
      name == "opp_expMayor" ~ "Mayor",
      name == "opp_expIncumbent" ~ "Incumbent",     
      name == "opp_ideolSomewhatideological" ~ "Somewhat ideological",
      name == "opp_ideolVeryideological" ~ "Very ideological"
    )) %>% 
  filter(name != "Intercept")

# add reference categories as 0
refcats2 <- tribble(~ attribute, ~ variable,
                   "Salary", "$15,000",
                   "Fundraising", "$25,000",
                   "Advertising", "No negative ads",
                   "Opponent's Experience", "No prior experience",
                   "Opponent's Ideology", "Moderate")
refcats2$value <- 0
imces_pes <- bind_rows(imces_pes, refcats2)

# put factors in the right order
imces_pes <- imces_pes %>% 
  mutate(attribute = factor(attribute, c("Salary", "Fundraising", "Advertising", 
                                         "Opponent's Experience", "Opponent's Ideology")),
         variable = factor(variable, c(
           "$15,000",
           "$50,000",
           "$80,000",
           "$25,000",
           "$100,000",
           "$300,000",
           "No negative ads",
           "Mostly negative ads",
           "No prior experience",
           "Mayor",
           "Incumbent",
           "Moderate",
           "Somewhat ideological",
           "Very ideological"
         )))
  


# make density plot. shows distribution of individual level MCEs
ggplot(imces_pes, aes(y = fct_rev(variable), x = value)) + 
  stat_halfeye(normalize = "groups", .width = c(.5, .9)) + 
  geom_vline(xintercept = 0, lty = 3) + 
  # geom_errorbarh(height = 0) + 
  # geom_point() + 
  facet_wrap(~attribute, ncol = 1, scales = "free_y") + 
  labs(x = "Individual-Level Effect", y = NULL) +
  ggtitle("Individual-Level Effects on\nInterest in Running") +
  scale_x_continuous(breaks = seq(-.4, .4, .1))
ggsave(file.path("figs/interest_conjoint_imce.pdf"), width=6*.75, height=8*.75)
ggsave(file.path(overleaf.dir, "figs/interest_conjoint_imce.pdf"), width=6*.75, height=8*.75)


summary_table <- imces_pes %>% 
  group_by(attribute, variable) %>% 
  summarise(AMCE = mean(value),
            `5th %` = quantile(value, .05),
            `25th %` = quantile(value, .25),
            `75th %` = quantile(value, .75),
            `95th %` = quantile(value, .95),
            .groups = "drop"
  ) 
summary_table %>% 
  kbl(
    format = "latex",
    booktabs = TRUE,
    digits = 3,
    linesep = "",
    caption = "Distribution of estimated IMCEs on interest in running.",
    label = "imce-interest"
  ) %>%
  group_rows(index = table(summary_table$attribute)) %>%
  kable_styling(latex_options = c("hold_position"), table.envir = FALSE) %>% 
  cat(file = file.path(overleaf.dir, "tables/imce-interest-summary.tex"))

# store results
imces_int <- imces_pes


# Results for win probability outcome --------------------------------------------


# extract point estimates of individual-level effects
imces_pes <- coef(winprob_reg_hier)[[1]][,1,] # point estimates of IMCEs
rids <- rownames(imces_pes) 
imces_pes <- as_tibble(imces_pes) %>% 
  mutate(responseid = rids) %>% 
  pivot_longer(-responseid) 

# tidy data frame for plotting  
imces_pes <- imces_pes %>% 
  mutate(
    attribute = case_when(
      grepl("salary", name) ~ "Salary",
      grepl("fundraising", name) ~ "Fundraising",
      grepl("opp_exp", name) ~ "Opponent's Experience",
      grepl("opp_ideol", name) ~ "Opponent's Ideology",
      grepl("advertising", name) ~ "Advertising"),
    variable = case_when(
      name == "Intercept" ~ "Intercept",
      name == "salary$50000" ~ "$50,000",
      name == "salary$80000" ~ "$80,000",
      name == "fundraising$100000" ~ "$100,000",        
      name == "fundraising$300000" ~ "$300,000",
      name == "advertisingMostlynegativeads" ~ "Mostly negative ads",
      name == "opp_expMayor" ~ "Mayor",
      name == "opp_expIncumbent" ~ "Incumbent",     
      name == "opp_ideolSomewhatideological" ~ "Somewhat ideological",
      name == "opp_ideolVeryideological" ~ "Very ideological"
    )) %>% 
  filter(name != "Intercept")

# add reference categories as 0
refcats2 <- tribble(~ attribute, ~ variable,
                    "Salary", "$15,000",
                    "Fundraising", "$25,000",
                    "Advertising", "No negative ads",
                    "Opponent's Experience", "No prior experience",
                    "Opponent's Ideology", "Moderate")
refcats2$value <- 0
imces_pes <- bind_rows(imces_pes, refcats2)

# put factors in the right order
imces_pes <- imces_pes %>% 
  mutate(attribute = factor(attribute, c("Salary", "Fundraising", "Advertising", 
                                         "Opponent's Experience", "Opponent's Ideology")),
         variable = factor(variable, c(
           "$15,000",
           "$50,000",
           "$80,000",
           "$25,000",
           "$100,000",
           "$300,000",
           "No negative ads",
           "Mostly negative ads",
           "No prior experience",
           "Mayor",
           "Incumbent",
           "Moderate",
           "Somewhat ideological",
           "Very ideological"
         )))



# make density plot. shows distribution of individual level MCEs
ggplot(imces_pes, aes(y = fct_rev(variable), x = value)) + 
  stat_halfeye(normalize = "groups", .width = c(.5, .9)) + 
  geom_vline(xintercept = 0, lty = 3) + 
  # geom_errorbarh(height = 0) + 
  # geom_point() + 
  facet_wrap(~attribute, ncol = 1, scales = "free_y") + 
  labs(x = "Individual-Level Effect", y = NULL) +
  ggtitle("Individual-Level Effects on\nWin Probability") +  
  scale_x_continuous(breaks = seq(-40, 40, 10))

ggsave(file.path(overleaf.dir, "figs/winprob_conjoint_imce.pdf"), width=6*.75, height=8*.75)
ggsave(file.path("figs/winprob_conjoint_imce.pdf"), width=6*.75, height=8*.75)

summary_table <- imces_pes %>% 
  group_by(attribute, variable) %>% 
  summarise(AMCE = mean(value),
            `5th %` = quantile(value, .05),
            `25th %` = quantile(value, .25),
            `75th %` = quantile(value, .75),
            `95th %` = quantile(value, .95),
            .groups = "drop"
  ) 
summary_table %>% 
  kbl(
    format = "latex",
    booktabs = TRUE,
    digits = 3,
    linesep = "",
    caption = "Distribution of estimated IMCEs on self-assessed probability of winning.",
    label = "imce-winprob"
  ) %>%
  group_rows(index = table(summary_table$attribute)) %>%
  kable_styling(latex_options = c("hold_position")) %>% 
  cat(file = file.path(overleaf.dir, "tables/imce-winprob-summary.tex"))


# store results
imces_winprob <- imces_pes





# compare indl-level effects ----------------------------------------------

comp <- imces_int %>% 
  select(responseid, name, value) %>% 
  filter(grepl("fundraising|opp_exp", name)) %>% 
  pivot_wider(names_from = "name",
              values_from = "value")
comp <- comp %>% 
  mutate(lofr_ratio = `fundraisingUSD100000` / opp_expIncumbent,
         hifr_ratio = `fundraisingUSD300000` / opp_expIncumbent)
quantile(comp$hifr_ratio, 1:20/20)
mean(comp$hifr_ratio >  1)

quantile(comp$lofr_ratio, 1:20/20)
mean(comp$lofr_ratio >  .5)



ggplot(comp) + 
  geom_density(aes(x = lofr_ratio, fill = "$100,000 fundraising"), alpha = .6) + 
  geom_density(aes(x = hifr_ratio, fill = "$300,000 fundraising"), alpha = .6) + 
  scale_x_continuous(limits = c(-.2, 3.5)) + 
  scale_fill_colorblind(name = NULL) + 
  labs(x = "Ratio of Fundraising IMCE to Incumbent IMCE") + 
  theme(legend.position.inside = c(.8, .9), legend.position = "inside")
ggsave(file.path(overleaf.dir, "figs/interest_imce_ratios.pdf"), width=6, height=4)
ggsave(file.path("figs/interest_imce_ratios.pdf"), width=8, height=6)

